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ABSTRACT 

In  this  paper,  we  compare  three  interpolation 
functions  in  a  discretized  continuum  when  used  in 
coupled  dynamic  atomistic-to-continuum  simulations. 
The  focus  is  on  assessing  the  ability  of  the  discrete 
continuum  model  to  capture  and  accurately  represent 
transient  effects,  namely  a  travelling  longitudinal  wave, 
through  both  the  mixed  atomistic-continuum  interface  and 
the  non-uniform  continuum  mesh  beyond.  We 
specifically  examine  the  differences  among  Bubnov- 
Galerkin,  partition  of  unity,  and  moving  least  squares 
finite  element  methods  in  the  continuum  part  of  the 
multiscale  model.  Our  study  shows  that  using  partition  of 
unity  interpolation  functions  in  the  continuum  produces 
superior  results  compared  to  the  other  two  approaches. 

INTRODUCTION 


In  this  paper  we  show  a  comparative  study  of  the 
quality  of  interpolation  that  best  suits  continuum  methods 
in  regions  at  and  near  the  interface  with  a  molecular 
dynamics  region.  We  specifically  examine  interpolation 
functions  prominent  in  general  finite  element  methods 
and  meshless  methods  -  Bubnov-Galerkin,  partition  of 
unity  [7],  and  moving  least  squares  [8]  -  and  assess  their 
ability  to  capture  a  travelling  wave  through  a 
discrete/continuum  interface  and  a  non-uniform  finite 
element  mesh  (increasing  element  size  away  from  the  MD 
region).  Within  the  interface  region,  where  the  continuum 
and  atomistic  scales  overlap,  the  displacements  on  the 
continuum  are  dictated  by  the  atomistic  results  generated 
from  MD.  In  this  study,  the  forces  between  the  domains 
are  communicated  from  the  atoms  to  the  continuum 
through  ghost  nodes. 

CONTINUUM  FORMULATION 


It  is  well  known  that  continuum  based  techniques 
such  as  Lagrangian  or  Eulerian  numerical  methods,  which 
use  constitutive  relations  that  do  not  account  for  the 
atomistic  structure,  are  invalid  beyond  the  scope  of  their 
calibration.  In  regions  containing  dislocations,  mobile 
defects,  or  nonlinear  material,  these  numerical  methods 
have  to  be  modified  to  capture  important  phenomena. 
Molecular  dynamics  (MD)  is  an  excellent  means  for 
predicting  interactions  on  an  atomic  scale  as  well  as 
predicting  the  response  when  sub-micron  scale 
phenomena  occur.  However,  MD  can  be  computationally 
expensive  beyond  small  sample  sizes  and  has  difficulty 
implementing  boundary  conditions  applied  at  a  continuum 
scale.  Therefore,  to  alleviate  these  problems  multiscale 
methods  have  been  developed  in  recent  years  to  couple 
the  continuum  and  atomistic  scales  together. 

There  has  been  extensive  work  on  developing  novel 
coupling  techniques  for  linking  atomistic  and  continuum 
scales.  These  techniques  include  the  quasicontinuum 
method  [1],  bridging  domain  method  [2],  bridging  scale 
method  [3]  and  homogenization  techniques  [4,5],  among 
others.  A  thorough  review  of  several  recent  techniques  is 
given  in  [6].  These  techniques  have  been  developed  using 
the  finite  element  method  within  the  continuum  scale. 
Though  seemingly  well  known,  to  our  knowledge,  an 
examination  of  the  level  of  approximation  and  choice  of 
interpolation  in  the  continuum  region  in  and  around  the 
discrete  atomistic  domain  has  not  been  shown. 


We  begin  by  reviewing  the  governing  equations  on 
the  continuum  scale.  The  conservation  of  momentum  can 
be  defined  as: 

Vo(LoP)+PoTofo-PoToii  (1) 

where  P  is  the  first  Piola-Kirchoff  stress  tensor,  fg  is  the 
body  force,  po  is  the  density,  Vg  is  the  initial  volume  and 
U  is  the  acceleration. 

From  classical  hyperelastic  continuum  approximation  we 
can  define  the  first  Piola-Kirchoff  stress  as: 


Fo  5F 


(2) 


where  W  is  the  potential  energy  density,  and  F  is  the 
deformation  gradient  defined  as: 


Sx  _  5u 
dX  ~  dX 


(3) 


where  X  denotes  the  reference  configuration  and  x 
denotes  the  spatial  or  current  configuration.  In  order  to 
use  equation  (1)  for  numerical  techniques  such  as  general 
finite  elements  we  use  the  principal  of  virtual  work  to 
obtain  the  variational  form: 


f w[v 0 (FoP)  +  PoVoh  -  =  0  (4) 

JQ.0 


where  w  is  the  virtual  displacement.  In  the  next  two 
sections  we  define  two  different  approaches  to 
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approximating  the  displacement  and  virtual  displacement, 
such  that: 


N 

u{x)  =  '^hi{x)ai  (5) 

/ 

where  h  is  a  vector  of  interpolation  functions  and  a  is  a 
vector  of  coefficients. 

PARTITION  OF  UNITY 

For  the  partition  of  unity  paradigm  [6],  define  a 
weighting  function,  W,  on  each  node  that  is  compactly 
supported  on  with  the  following  properties: 

1.  JVj(x)eC^(Bj}  s>0_ 

2.  JVj(x)>0  VxeBj 

3.  Wi{x)  =  0  elsewhere 

The  symbol  Cq  [Bj  )  stands  for  the  space  of  functions  that 
are  compactly  supported  on  Bj ,  where  in  the  case  of 
general  finite  elements  Bj  is  generated  using  neighboring 

elements,  which  have  continuous  derivatives  of  order  s. 
We  define  the  Shepard  partition  of  unity  function  at  each 
node  /  as: 


^/(x) 


(6) 


From  the  partition  of  unity  property  it  follows  that  the 
functions  satisfy  zeroth  order  consistency,  i.e.  they  ensure 
that  rigid  body  modes  are  exactly  satisfied.  The  next  step 
is  to  develop,  at  each  node  I,  a  local  approximation  space 

=  span^^^  {p^  (x)}  c;  // '  (5^  n  n)  (7) 


where  h  is  a  measure  of  the  size  of  the  spheres,  p  is  the 
polynomial  order,  ^  is  an  index  set,  H'  is  the  first  order 
Hilbert  space,  and  /?^(x)  is  a  polynomial  or  other 
function.  Finally,  the  global  approximation  space  is 
defined  by  pasting  together  the  local  spaces  as  follows: 

V’'’”  (8) 


Hence,  any  function  e  can  be  written  as: 

^  X  S  *1” 

/=1 

hM  =  (P^!{^)pA^)  (10) 

and  him  is  the  shape  function  at  node  1  corresponding  to 
the  m*  degree  of  freedom. 


In  moving  least  squares  we  set  the  approximation  to: 

(x)  =^Pm  (x)  =  (x)p(x)  (11) 

m=l 

where  p  is  a  vector  composed  of  the  monomial  basis 
functions  as  in  equation  (7)  and  P(x)  is  a  vector  composed 
of  their  coefficients.  These  coefficients  are  obtained  by 
using  a  weighted  least  square  fit  for  the  local 
approximation.  We  can  derive  this  by  minimizing  the 
difference  between  the  local  approximation  and  the 
function,  such  that: 


|^  =  A(x)p(x)-B(x)a  =  0 

(12) 

a(x)=P^T'(x)p 

(13) 

B(x)  =  P^T'(x) 

(14) 

where  P  is  a  matrix  composed  of  the  monomial  basis 
functions  and  ^  is  a  matrix  composed  of  the  weighting 
functions  having  the  same  properties  as  those  used  in 

partition  of  unity  interpolation  functions. 

This  results  in: 

v^’^(x)  =  p^(x)A  '(x)B(x)a  =  ^/i/(x)«^  (15) 

/ 


And  the  shape  function  is  defined  as: 

^/(x)  =  P^(x:)a^'(x)p(x;)1F;(x)  (16) 

MOLECULAR  DYNAMICS 

For  the  atomistic  scale  the  governing  equation  for 
MD  is  Newton’s  equation  of  motion  defined  as 

ma»a=fa  (17) 

where  nia  is  the  mass,  is  the  acceleration  anAfa  is  the 
force  acting  on  discrete  atoms,  a.  For  our  study  we  will 
only  examine  short  range  interactions.  The  force  is 
defined  as: 


ab 


b^a 


dx„ 


where  is  the  interatomic  potential. 


(18) 


In  this  paper  we  specifically  examine  a  linear 
harmonic  potential  and  a  non-linear  Lennard-Jones 
potential,  where  the  harmonic  potential  is  given  as: 

^l=\k{r^b-rabfif  (19) 


MOVING  LEAST  SQUARES 


where  A:  is  a  eonstant  and  ro  is  the  zero  potential  distanee 
between  two  atoms.  The  Lennard-Jones  potential  is 
defined  as: 


Vab 


f 

a 

Jab, 


where  s  and  o  and  constants. 


(20) 


performed  using  MD  throughout  the  entire  model  using 
79  atoms.  For  the  comparative  study,  the  domain  is 
discretized  such  that  for  -2  <  x  <  2  each  atom  is 
individually  resolved  and  from  -10  <  x  <  2  and 
2<x<10  the  different  numerical  interpolation  schemes 
are  applied.  We  compare  the  use  of  Bubnov-Galerkin, 
finite  elements,  partition  of  unity  and  moving  least 
squares  interpolation  functions  to  full  MD  throughout  the 
domain.  The  initial  wave  function  is  shown  in  figure  2. 


COUPLING 

In  figure  1  we  define  the  domain  as  discretized  into  a 
region  in  which  the  continuum  equations  are  applied,  O^, 
and  a  region  in  which  MD  is  applied,  Q^.  There  is  an 
overlap  between  these  two  regions  defined  as  the  interface 
region  Q'. 


Figure  1:  Coupling  Domain 

The  constraint  that  matches  atoms  to  nodes  in  the 
interface  is  applied  through  a  penalty  formulation.  The 
result  is  a  modified  variational  form  of  (4): 


^M  +  f;y(h(xffh(xf} 


/=1 

-f  i„t 

^  1  ^ 

—  M 


a  = 


(21) 


AC 


(2'a-'-'''a)+^/h(x/fu; 


i=i 


Where  is  the  penalty  constant  which  is  generally  a  large 
positive  number.  To  enforce  the  interface  conditions  for 
MD,  ghost  atoms  are  placed  in  the  continuum  region  (see 
figure  1)  to  avoid  a  surface  layer  or  unphysical 
termination  in  the  interface. 

NUMERICAL  EXAMPLES 

We  present  three  examples  of  a  Gaussian  wave 
propagating  through  atomic  medium  to  illustrate  the 
preliminary  results.  The  first  two  examples  are 
demonstrated  in  a  ID  domain,  while  the  third  is  example 
is  in  two  dimensions. 

For  the  first  two  examples  the  analysis  is  preformed 
on  a  ID  section  of  atoms  and  compared  with  simulation 


Figure  2:  ID  Example  formulation 

For  the  first  example  we  use  a  harmonic  interatomic 
potential.  The  results  of  the  error  in  the  displacements  as 
the  wave  propagates  are  shown  in  figure  3.  The  red  line 
represents  the  error  in  a  continuum  disctretized  using 
finite  elements  the  green  line  represents  moving  least 
squares  functions  and  the  blue  represents  partition  of 
unity  functions. 


Figure  3:  Error  comparisons  for  Harmonic  potential. 

From  this  result,  the  partition  of  unity  interpolation 
scheme  is  perform  extremely  well  when  compared  to  the 
fiill-MD  analysis.  Whereas  finite  elements  and  MLS 
functions  perform  about  the  same.  Spikes  in  the  error 


occur  for  wall  three  techniques  when  the  wave  passes 
through  the  boundaries. 

For  the  second  example  we  compare  the  analysis 
shown  in  figure  2  using  a  Lennard- Jones  potential.  The 
results  are  shown  in  figure  4.  The  results  colaberate  the 
above  example  with  partition  of  unity  functions 
performing  better  then  the  other  two  techniques. 


Figure  4:  Error  comparisons  for  Lennard-Jones 
potential 

For  the  2D  example  the  domain  is  discretized  such 
that  for  -2<  x,y  <2  each  atom  is  individually  resolved 
and  from  -10<x,y<2  and  2  <  x,  y  <  10  the  different 
numerical  interpolation  schemes  are  applied.  We  compare 
the  use  of  Bubnov-Galerkin,  finite  elements  and  partition 
of  unity  interpolation  functions  to  full  MD  throughout  the 
domain.  We  use  a  harmonic  interatomic  potential  in  this 
example.  The  initial  wave  function  is  shown  in  figure  5. 


Figure  5:  2D  example  formulation 

The  results,  shown  in  figure  6,  demonstrate  again  that 
partition  of  unity  shape  function  out  perform  finite 
elements,  however,  the  margin  is  much  less. 


Figure  3:  Error  comparisons  for  2D  example. 
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